library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ ggplot2 3.4.0     ✔ purrr   0.3.4
## ✔ tibble  3.1.7     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(here)
## here() starts at /Users/jasonmoggridge/Documents/CUBE_UOttawa
library(ggiraph)
library(patchwork)

source(here('R/report_plots.R'))
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
swabs <- read_rds(here('data/cube.rds')) |> 
  filter(str_detect(site, '^UO_')) |> 
  glimpse()
## Rows: 721
## Columns: 39
## $ date_swab              <date> 2021-09-21, 2021-09-21, 2021-09-21, 2021-09-21…
## $ city                   <fct> Ottawa, Ottawa, Ottawa, Ottawa, Ottawa, Ottawa,…
## $ site                   <fct> UO_90U, UO_90U, UO_90U, UO_DMS, UO_DMS, UO_DMS,…
## $ floor                  <chr> "Level 1", "Level 1", "Level 1", "Level 1", "Le…
## $ location               <chr> "Site 1: Entrance Foyer, Just Before The Securi…
## $ replicate              <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ lot                    <chr> "XG024", "XG024", "XG024", "XG024", "XG024", "X…
## $ barcode                <chr> "10349", "2921", "10310", "14573", "11841", "10…
## $ pcr_positive           <fct> Negative, Negative, Negative, Negative, Negativ…
## $ worker                 <fct> Aaron, Alex, Rees, Partha, Naila, Keaton, Rache…
## $ lab                    <fct> Ottawa, Ottawa, Ottawa, Ottawa, Ottawa, Ottawa,…
## $ type                   <fct> University, University, University, University,…
## $ outbreak               <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ susp_outbreak          <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ outbreak_flag          <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ comments               <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ room_demographics      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ cases                  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ case_dates             <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ case_comment           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ date_pcr               <date> 2021-09-24, 2021-09-24, 2021-09-24, 2021-09-24…
## $ pcr_result             <chr> "not detected", "not detected", "not detected",…
## $ negative_control       <lgl> FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, …
## $ spoiled                <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
## $ co2                    <dbl> 538, 519, 538, 582, 528, 582, 730, 687, 730, 54…
## $ date_shipped           <date> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ pcr_result2            <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ pcr_result3            <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ location_within_room   <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ comment                <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ case_comments          <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ redo                   <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ swab_type              <chr> "flock", "flock", "flock", "flock", "flock", "f…
## $ case_role              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ pcr_ct                 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ floor_location         <chr> "Level 1: Site 1: Entrance Foyer, Just Before T…
## $ sample_location        <fct> "Uo_90u: Level 1: Site 1: Entrance Foyer, Just …
## $ replicate_barcode      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ replicate_swab_1_dates <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
swabs |> count(site, date_swab)
## # A tibble: 316 × 3
##    site   date_swab      n
##    <fct>  <date>     <int>
##  1 UO_90U 2021-09-21     3
##  2 UO_90U 2021-09-23     3
##  3 UO_90U 2021-09-28     3
##  4 UO_90U 2021-09-30     3
##  5 UO_90U 2021-10-05     2
##  6 UO_90U 2021-10-07     2
##  7 UO_90U 2021-10-12     2
##  8 UO_90U 2021-10-14     2
##  9 UO_90U 2021-10-19     2
## 10 UO_90U 2021-10-21     3
## # … with 306 more rows
# waste water data 
raw_ww <-
  readxl::read_excel(here::here('data/wastewater.xlsx')) |>
  janitor::clean_names()

uottawa_ww <-
  raw_ww |>
  filter(!is.na(sample_date)) |>
  mutate(
    source =
      source_name |>
      str_remove('^Data - uOttawa - ') |>
      str_remove('.xlsx$'),
    site = case_when(
      source == 'uo_aa' ~ 'AA',
      source == 'uo_fa' ~ 'FA',
      source == 'uo_ft' ~ 'FT',
      source == 'uo_na' ~ 'NA',
      source == 'uo_rgn' ~ 'RGN',
      source == 'uo_st' ~ 'ST',
      source == 'uo_tbt' ~ 'TBT',
      TRUE ~ 'not matched'
    )
  )

# ottawa wastewater data: summary
summarise_wastewater <- function(){
  read_rds(here('data/ww_ottawa.rds')) |> 
  select(sample_date, starts_with('cov_')) |> 
  pivot_longer(contains('cov_')) |> 
  mutate(target = str_extract(name, 'cov_n.'),
         stat = str_extract(name, 'mean|sd'),
         ) |> 
  select(-name) |> 
  pivot_wider(names_from = stat, values_from = value) |> 
  mutate(.lower = mean - sd, .upper = mean + sd)   
}

wifi <- 
  read_rds(here('data/uo_wifi.rds')) |> 
  filter(date >= min(swabs$date_swab), date <= max(swabs$date_swab))

📊 Summary

This plot shows the counts of positive (red) and negative (yellow) samples collected at each facility over time (x-axis). Samples that could not be tested are shown in navy. Only flocked swabs are shown. (Other sponge swabs were collected on 2022-04-28 were for comparing flocked and sponge swabs.)

plot_timeseries_summary(
  swabs |> 
    filter(
      swab_type %in% 'flock',
      !negative_control,
      !is.na(negative_control)
    ),
  height_svg = 3.3,
  width_svg = 8
  )

⡯⡷ Dotplot

This plot shows the counts of positive and negative samples collected at each facility over time.

plot_timeseries_locations(
  swabs |> 
    filter(swab_type == 'flock'),
  height_svg = 6,
  width_svg = 8
  )

Waste Water Data

uottawa_ww |>
  ggplot(aes(sample_date, viral_copies_per_copies_pep_avg, color = site)) +
  geom_line(alpha = 0.5) +
  geom_point(alpha = 0.5) +
  labs(x = 'Date', y = 'SARS-CoV-2 / PPMoV')
## Warning: Removed 16 rows containing missing values (`geom_line()`).
## Warning: Removed 80 rows containing missing values (`geom_point()`).

rm(raw_ww, uottawa_ww)

UO Results

This section contains results from modeling covid cases at UO using swab results as a predictor.


Modelling

Derek and Jason used mixed effects logistic regression to model the occurrence of cases (y/n) based on swab results for the previous week (the proportion of positive swabs). The model has a random intercept for sites.

Model specification

I switched this model to a Bayesian mixed-effects model. In the mixed-effects generalized linear model, MRT/90U are problem-sites where the random intercept parameter causes the model fit to be overspecified (model matrix becomes singular).

Our model formula is cases ~ positives[lag 1 week] + (1|site).

The model is fit as follows:

swab_model <-
  blme::bglmer(
    cases_binary ~ detection_lag_1week + (1 | site),
    data =  uo_sites,
    family = binomial
  )
swab_model <- read_rds(here('model/uo_swab_model.rds'))


Model response time-series

These plots show the swab results, cases, and predictions by the current model.

# uo cases data - long table
cases_long <-
  read_rds(here('./data/latest_uo_cases.rds')) |>
  filter(!is.na(positive_result),
         positive_result > lubridate::as_date('2021-09-14')) |>
  select(case, role, confirmed_by_test, positive_result, c(`90U`:TBT)) |>
  mutate(
    `Campus Total` = TRUE,
    # 3 days transmissibility
    transmission_start = positive_result - days(3)
    ) |>
  pivot_longer(`90U`:`Campus Total`, names_to = 'site') |>
  filter(value) |>
  select(-value) |> 
  mutate(site = fct_relevel(site, "Campus Total", after = 0L))

uo_pred <- read_rds(here('model/uo_pred.rds'))

# create a 3-panel fig for each of 6 buildings.
map(
  # site names for filtering 
  .x = swabs$site |> unique() |> str_remove('UO_'),
  .f = ~plot_uo_pred_swab_case(
    swabs = swabs,
    cases = cases_long,
    preds = uo_pred,
    site_select = .x
  )
) |> wrap_plots(guides = 'collect', ncol = 2)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.


Model summary

These tables show the model coefficients and statistics.

broom.mixed::glance(swab_model) |> 
  mutate_if(is.numeric, round, 2) |> 
  kableExtra::kable(caption = 'Model statistics',
                    format = "html") |> 
  kableExtra::kable_styling(bootstrap_options = "striped", 
                            full_width = F, position = "left")
## Loading required package: blme
## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
Model statistics
nobs sigma logLik AIC BIC deviance df.residual
84 1 -24.52 55.04 62.33 42.83 81 
broom.mixed::tidy(swab_model) |>
  mutate(p.value = as.character(p.value)) |> 
  mutate_if(is.numeric, round, 3) |> 
  mutate(p.value = as.numeric(p.value)) |> 
  kableExtra::kable(caption = 'Model parameters',     
                    format = "html") |> 
  kableExtra::kable_styling(bootstrap_options = "striped",
                            full_width = F, position = "left")
Model parameters
effect group term estimate std.error statistic p.value
fixed NA (Intercept) -3.301 0.804 -4.107 0.0000401
fixed NA detection_lag1w 7.114 2.112 3.369 0.0007549
ran_pars site sd__(Intercept) 1.067 NA NA NA 
# tidy(swab_model, effects = 'ran_vals')


sjPlot::tab_model(
  swab_model,
  show.re.var = TRUE,
  pred.labels = c("(Intercept)", "Swabs @ t-1week"),
  dv.labels = "Relation between UO cases (y/n) and proportion\n of positive swabs the previous week"
)
  Relation between UO cases (y/n) and proportion of positive swabs the previous week
Predictors Odds Ratios CI p
(Intercept) 0.04 0.01 – 0.18 <0.001
Swabs @ t-1week 1229.31 19.59 – 77128.28 0.001
Random Effects
σ2 3.29
τ00 site 1.14
ICC 0.26
N site 6
Observations 84
Marginal R2 / Conditional R2 0.378 / 0.538


Model response curves

This plot shows the modelling data as points (detection lagged 1 week and cases at a given site- y/n), as well as how the probability of future cases varies by the previous weeks detection level and site, according to our model (curves).

read_rds(here('fig/plt_model_resp_curves.rds')) +
  labs(subtitle = 'Points show case/swab data and curves show model probability values')
## Warning: Removed 42 rows containing missing values (`geom_point()`).


Random Effects for Sites

This plot shows the odds ratios for the intercepts of the model (an intercept for each site).

# random effect plot: site intercepts; axes are flipped
sjPlot::plot_model(swab_model, type = 're') +
  theme_light() +
  labs(title = "Random effects", y = 'Odds Ratio')


Swabs and Cases

This panel shows the cases counts at UO over the course of our sampling period. The case data shown represents the days on which a positive test was reported (black rug lines) and the presumed start of transmissiblity for each case (red lines).

# generated from file R/plot_UO_case_load.R
read_rds(here('fig/plt_cases_swabs.rds'))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 4 rows containing missing values (`geom_bar()`).


Swabs, Wifi, and CO2

This panel shows linked data from swab results, CO2 readings, and wifi traffic (number of users @ peak, daily) during our study period. Unfortunately, we do not have wifi data for 90U. 

source(here('R/uott_presentation.R'))
uo_multivariate & theme(text = element_text(size = 13))
## Warning: Removed 1 rows containing missing values (`geom_point()`).

rm(uo_multivariate)


Wifi Traffic

plot_UO_wifi_ts(wifi, swabs)

This panel shows a time-series of the daily peak number of wifi users at UO facilities. Sampling days are highlighted in blue.